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Abstract. We implement a two-stage approach of the Wang-Landau algorithm to investigate the critical 
properties of the 3D Ising model with quenched bond randomness. In particular, we consider the case where 
disorder couples to the nearest-neighbor ferromagnetic interaction, in terms of a bimodal distribution of 
strong versus weak bonds. Our simulations are carried out for large ensembles of disorder realizations and 
lattices with linear sizes L in the range L = 8 — 64. We apply well-established finite-size scaling techniques 
and concepts from the scaling theory of disordered systems to describe the nature of the phase transition 
of the disordered model, departing gradually from the fixed point of the pure system. Our analysis (based 
on the determination of the critical exponents) shows that the 3D random-bond Ising model belongs to 
the same universality class with the site- and bond-dilution models, providing a single universality class 
for the 3D Ising model with these three types of quenched uncorrelated disorder. 



PACS. PACS. 05.50-|-q Lattice theory and statistics (Ising, Potts, etc. 
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1 Introduction 

Understanding the role of impurities on the nature of 
phase transitions is of great importance, both from ex- 
perimental and theoretical point of view W. First-order 
phase transitions are known to be significantly softened 
under the presence of quenched randomness [2 3 4 5 6 7 8, 
whereas continuous transitions may have their exponents 
altered under random fields or random bonds |9I10) . There 
are some very useful phenomenological arguments and 
some, perturbative in nature, theoretical results, pertain- 
ing to the occurrence and nature of phase transitions un- 
der the presence of quenched randomness [4I5I11I12] . His- 
torically, the most celebrated criterion is that suggested by 
Harris [5]. This criterion relates directly the persistence, 
under random bonds, of the non random behavior to the 
specific heat exponent ap of the pure system. According 
to this criterion, if Up > 0, then disorder will be rele- 
vant, i.e., under the effect of the disorder, the system will 
reach a new critical behavior. Otherwise, if ap < 0, disor- 
der is irrelevant and the critical behavior will not change. 
Pure systems with a zero specific heat exponent (ap = 0) 
are marginal cases of the Harris criterion and their study, 
upon the introduction of disorder, has been of particular 
interest [13] . The paradigmatic model of the marginal case 
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is, of course, the general random 2D Ising model, which 
has been extensively debated [13]. 

Respectively, the 3D Ising model with quenched ran- 
domness - which is a clear case in terms of the Harris cri- 
terion having a positive specific heat exponent in its pure 
version - has also been extensively studied using Monte 
Carlo (MC) simulations 15 16 17 18 19 20 21 22 23 24j|2S] 
and field theoretical renormalization group approaches [26127128] . 
Especially, the diluted model can be treated in the low- 
dilution regime by analytical perturbative renormaliza- 
tion group methods [29'30'31!, where a new fixed point 
independent of the dilution has been found, yet for the 
strong dilution regime only MC results remain valid. Al- 
though the first numerical studies of the model suggested 
a continuous variation of the critical exponents along the 
critical line, it soon became clear that the concentration- 
dependent critical exponents found in MC simulations are 
the effective ones characterizing the approach to the asymp- 
totic regime |17I18I19] . Note, here, that a crucial problem 
of the new critical exponents obtained in these studies is 
that the ratios P/v and -f/i' occurring in finite-size scal- 
ing (FSS) analysis are almost identical for the disordered 
and pure models. In fact, for the pure 3D Ising model, ac- 
curate values are [32]: = 0.6304(13), P/i' = 0.517(3), 
7/1/ = 1.966(3), and a = 0.1103(1). Respectively, for 
the site- and bond-diluted models, the most accurate sets 
of asymptotic exponents {v, P/v, and 7/1^) have been 
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given by the numerical works of Ballesteros et al. ^I] ^'^d 
Berche et al. ^ are: (0.6837(53), 0.519(3), 1.963(5)) and 
(0.68(2), 0.515(5), 1.97(2)). 

The above estimates of critical exponents provided ev- 
idence that the 3D Ising model with quenched uncorre- 
cted disorder belongs to a single universality class, dis- 
tinct from that of the pure model, as also indicated by 
the Harris criterion, independent of the considered disor- 
der distribution. Yet, in a very recent paper Murtazaev 
and Babaev [33] using MC simulations and FSS meth- 
ods on the site-diluted model, the above view was con- 
tradicted and these authors suggested that the model has 
two regimes of critical behavior universality, depending 
on the nonmagnetic impurity concentration. Motivated by 
the above contradictions and the great theoretical interest 
of the existence of universality classes in disordered mod- 
els, it seems favorable to investigate the 3D Ising model 
with bond disorder in order to compare all these three 
kinds of disorder (site-, bond-dilution, and bond disorder) 
and to verify whether these lead to the same set of new 
critical exponents, as would be, in principle, expected by 
universality arguments [35]. To this end, the first approach 
has been recorded in a recent brief report by the present 
authors [34, who studied the random-bond model for a 
single value of the disorder strength, obtaining values for 
the critical exponents in very good agreement with the 
most accurate estimates for the bond- and site-diluted 
models. Here, we extend the analysis of reference [M] , 
including more values of the disorder strength and sev- 
eral aspects of the FSS behavior of the model, thus pre- 
senting a complete picture of the disorder-induced phase 
transition of the 3D random-bond Ising model. The main 
outcome of our work is that, indeed, the 3D Ising model 
with quenched, uncorrelated bond disorder belongs to the 
same universality class as the site- and bond-diluted mod- 
els, defining in this way a complete universality class in 
disordered spin models. 

The rest of the paper is organized as follows: In Sec- 
tion [2] we define the model by implementing the bond- 
disorder distribution and describe the basic elements of 
our two-stage numerical approach. In Section[3]we present 
a detailed FSS analysis of the obtained numerical data, es- 
timate critical exponents with high accuracy, and discuss 
the disorder-induced second-order phase transition of the 
model under the general prism of universality. This con- 
tribution closes with Section |4l where a brief summary 
of our conclusions and an outlook for further studies are 
given. 



2 Model and simulation method 

We start this Section with the definition of the model. 
In the following we consider the 3D bond-disorder Ising 
model, whose Hamiltonian with uncorrelated quenched 
random interactions reads 

H = -^J^jSiSj, (1) 



where the spin variables Si take on the values — 1,+1, 
(ij) indicates summation over all nearest-neighbor pairs 
of sites, and the ferromagnetic interactions Jij > follow 
a bimodal distribution of the form 

P{J^,) - \ [SiJ^J - Jl) + 5{Jn - J2)]- (2) 

In equation ^ we set J1 + J2 = 2, Ji > J2, and r = J2/J1 
reflects the strength of the bond randomness. Addition- 
ally, in the following, we fix as usual 2/cb/(Ji + J2) — 1, 
to set the temperature scale (ks = 1 also for simplicity). 
The values of the disorder strength r considered through- 
out this paper are the following: r = 0.75/1.25, 0.5/1.5, 
and 0.25/1.75. 

Resorting to large scale MC simulations is often nec- 
essary, especially for the study of the critical behavior of 
disordered systems. It is also well known that for such 
complex systems traditional methods become inefficient 
and thus in the last few years several sophisticated al- 
gorithms, some of them are based on entropic iterative 
schemes, have been proven to be very effective [HSj. In the 
last few years we have used an entropic sampling imple- 
mentation of the Wang-Landau (WL) algorithm |36l37j 
to study some simple [38139] . but also some more com- 
plex systems [40141142143144] . One basic ingredient of this 
implementation is a suitable restriction of the energy sub- 
space for the implementation of the WL algorithm. This 
was originally termed as the critical minimum energy sub- 
space restriction ,38 39 and it can be carried out in many 
alternative ways, the simplest being that of observing the 
finite-size behavior of the tails of the energy probability 
density function of the system [35] . 

Complications that may arise in complex systems, i.e., 
random systems or systems showing a first-order phase 
transition, can be easily accounted for by various simple 
modifications that take into account possible oscillations 
in the energy probability density function and expected 
sample-to-sample fluctuations of individual realizations. 
Recently, details of various sophisticated routes for the 
identification of the appropriate energy subspace {Ei, E2) 
for the entropic sampling of each realization have been 
presented in references [41142] . To estimate the appropri- 
ate subspace from a chosen pseudocritical temperature 
one should be careful to account for the shift behavior 
of other important pseudocritical temperatures and ex- 
tend the subspace appropriately from both low- and high- 
energy sides (as also discussed in reference [H]) in order 
to achieve an accurate estimation of all finite-size anoma- 
lies. Of course, taking the union of the corresponding sub- 
spaces, ensures accuracy for the temperature region of all 
studied pseudocritical temperatures. 

The up to date version of our implementation uses a 
combination of several stages of the WL process. First, we 
carry out a starting (or preliminary) multi-range (multi- 
R) stage, in a very wide energy subspace. This prelimi- 
nary stage is performed up to a certain level of the WL 
random walk. The WL refinement is G{E) — ^ / • G{E), 
where G{E) is the density of states (DOS) and we fol- 
low the usual modification factor adjustment fj+i = y/Jj 
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Fig. 1. Disorder distribution of the susceptibility maxima of a 
lattice with linear size L = 32 for disorder strength r = 0.5/1.5. 
The running average over the samples is shown by the thick 
solid line. 

and /i = e [38139) . The preliminary stage may consist of 
the levels : j = = 18 and to improve accuracy 

the process may be repeated several times. However, in 
repeating the preliminary process and in order to be effi- 
cient, we use only the levels j — 13, . . . , 18 after the first 
attempt, using as starting DOS the one obtained in the 
first random walk at the level j = 12. From our experi- 
ence, this practice is almost equivalent to simulating the 
same number of independent WL random walks. Also in 
our recent studies we have found out that is much more 
efficient and accurate to loosen up the originally applied 
very strict flatness criteria [3S]. Thus, a variable flatness 
process starting at the first levels with a very loose flat- 
ness criteria and assuming at the level j = 18 the original 
strict flatness criteria is nowadays used. After the above 
described preliminary multi-R stage, in the wide energy 
subspace, one can proceed in a safe identification of the 
appropriate energy subspace using one or more alterna- 
tives outlined in reference [391 . 

The process continues in two further stages (two-stage 
process), using now mainly high iteration levels, where 
the modification factor is very close to unity and there is 
not any significant violation of the detailed balance con- 
dition during the WL process. These two stages are suit- 
able for the accumulation of histogram data (for instance, 
energy- magnetization histograms), which can be used for 
an accurate entropic calculation of non-thermal thermo- 
dynamic parameters, such as the order parameter and its 
susceptibility [31]. In the first (high-level) stage, we follow 
again a repeated several times (typically ~ 5 — 10) multi- 
R WL approach, carried out now only in the restricted 
energy subspace. The WL levels may be now chosen as 
j = 18, 19, 20 and as an appropriate starting DOS for the 
corresponding starting level the average DOS of the pre- 
liminary stage at the starting level may be used. Finally, 
the second (high-level) stage is applied in the refinement 
WL levels j = ji, . . . ,ji + 3 (typically ji = 21), where 
we usually test both an one-range (one-R) or a multi- 



Fig. 2. (Color online) Disorder-averaged specific heat [C]av 
(lower curves) and magnetic susceptibility [x]av (upper curves) 
as a function of the temperature T for a lattice with linear 
size L = 8. The ensemble of random realizations is Q = 2000 
for all cases of the disorder strength shown here. Notice the 
expected shift of the pseudocritical temperatures to the left 
with increasing the disorder strength (r — >■ 0). 

R approach with large energy intervals. In the case of 
the one-R approach we have fotmd very convenient and 
more accurate to follow the Belardinelli-Pereyra [4Bj ad- 
justment of the WL modification factor according to the 
rule In/ ^ t~^. Finally, it should be also noted that by 
applying in our scheme a separate accumulation of his- 
togram data in the starting multi-R stage (in the wide 
energy subspace) offers the opportunity to inspect the be- 
havior of all basic thermodynamic functions in an also 
wide temperature range and not only in the neighborhood 
of the finite-size anomalies. 

Using this scheme we performed extensive simulations 
for several lattice sizes in the range L ~ 8 — 64, over 
large ensembles {1, ■••,(/,•••, Q} of random realizations 
{Q = 1000 — 3000). It is well known that, extensive dis- 
order averaging is necessary for the study of random sys- 
tems, where usually broad distributions are expected lead- 
ing to a strong violation of self-averaging [22146] . A mea- 
sure from the scaling theory of disordered systems, whose 
limiting behavior is directly related to the issue of self- 
averaging [22146] may be defined with the help of the rel- 
ative variance of the sample-to-sample fluctuations of any 
relevant singular extensive thermodynamic property Z as 
follows: Rz = {[ZXv - Figure [Jj presents 

evidence that the above number of random realizations is 
sufficient in order to obtain the true average behavior and 
not a typical one. In particular, we plot in this figure (for 
a lattice size L = 32 and disorder strength r = 0.5/1.5) 
the disorder distribution of the susceptibility maxima x* 
(where the subscript q denotes the random realization) 
and the corresponding running average, i.e., a series of 
averages of different subsets of the full data set - each of 
which is the average of the corresponding subset of a larger 
set of data points, over the samples for the simulated en- 
semble of Q = 1447 disorder realizations. A first striking 
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Fig. 3. Shift behavior of several pseudocritical temperatures 
defined in the text for r — 0.75/1.25. The error bars reflect the 
sample-to-sample fluctuations. 



Fig. 4. Same as Figure|3] but for r = 0.5/1.5. The inset shows 
the FSS of the sample-to-sample fluctuations of the pseudo- 
critical temperature of the magnetic susceptibility. 



observation from this figure is the existence of very large 
variance of the values of x*, indicating the violation of 
self-averaging for this quantity |22.46. . This figure illus- 
trates that the simulated number of random realizations 
is sufiicient in order to probe correctly the average behav- 
ior of the system, since already for Q « 300 the average 
value of Xq is stable. 

Closely related to the above issue of self-averaging in 
disordered systems is the manner of averaging over the dis- 
order [47] . This non-trivial process may be performed in 
two distinct ways when identifying the finite-size anoma- 
lies, such as the peaks of the magnetic susceptibility. The 
first way corresponds to the average over disorder real- 
izations ([. . .]av) and then taking the maxima ([. . .]*^,), or 
taking the maxima in each individual realization first, and 
then taking the average ([. . ■*]av)- In tbe present paper we 
present our FSS analysis using mainly the first approach 
of averaging, although we should note that also the sec- 
ond gave comparable results. As an example, we show in 
Figured] the curves of the specific heat and magnetic sus- 
ceptibility for a lattice with linear size L = 8 averaged over 
Q = 2000 random realizations. One may retrieve the loca- 
tion of the pseudocritical point by taking the maximum of 
these curves. Commenting on the statistical errors of our 
WL scheme, we found that the statistical errors of our 
scheme on the observed average behavior were of small 
magnitude (of the order of the symbol sizes) and thus are 
neglected in the figures. On the other hand, for the case 
[. . .*]av the error bars shown reflect the sample-to-sample 
fluctuations. 



3 Numerical results and discussion 

We present in this Section our numerical results and FSS 
analysis for the 3D random-bond Ising model. In FiguresjS] 
- |S] we illustrate in the main panels the shift behavior 
of several pseudocritical temperatures, i.e., we take the 



average over the pseudocritical temperatures of the in- 
dividual samples. In all cases the error bars reflect the 
sample-to-sample fluctuations. The considered pseudocrit- 
ical temperatures correspond to the peaks of the following 
six quantities: speciflc heat C, magnetic susceptibility x, 
derivative of the absolute order parameter with respect to 
inverse temperature K = 1/T [3B] 

d{\M\)/dK={\M\H)-m){H) (3) 

and logarithmic derivatives of the first {n = 1), second 
(n — 2), and fourth (n = 4) powers of the order parameter 
with respect to inverse temperature HI] 

dln{M")/dK = {M"H)/{M") - (H). (4) 

Fitting our data for the whole lattice range to the expected 
power-law behavior 

[T*zU=T, + bL-'/\ (5) 

where Z stands for the different thermodynamic quantities 
mentioned above, we estimate the critical temperatures 
as a function of the disorder strength and the respective 
critical exponents of the correlation length. In particular, 
for the case of very weak disorder (case r — 0.75/1.25 
shown in Figure [3]) we estimate a value Tc = 4.4873(21) 
for the critical temperature and ly = 0.6511(33) for the 
respective critical exponent. Although this value of v is 
larger than the corresponding value of the pure system, 
it is still far away from the expected value of the random 
system (as obtained in references 121125134] and Figures H] 
and |S]). This is due to finite-size effects that dominate 
in the regime of very weak randomness, hindering the ap- 
proach of the system to the asymptotic limit. On the other 
hand, the simultaneous fittings of the form ([5]) shown in 
the main panels of Figures U] and |S] for stronger values 
of the disorder strength r reveal estimates for the corre- 
lation length exponent around the value 0.6849(50), in 
excellent agreement with the values 0.6837(53) and 0.68(2) 
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Fig. 5. Same as Figure |4]but for the case r = 0.25/1.75. 

given by Ballesteros et al. fTT and Berche et al. [25] for 
the site- and bond-diluted models. Additionally, the esti- 
mated values for the critical temperatures, in all Figures O 
m and [5l follow the usual reduction trend with increasing 
disorder strength (r — >■ 0), as already shown in Figure [2] 
of Section [21 

Using the above sample-to-sample fluctuations of the 
pseudocritical temperatures and the theory of FSS in dis- 
ordered systems as introduced by Aharony and Harris [46| 
and Wiseman and Domany [22^, one may further exam- 
ine the nature of the fixed point controlling the critical 
behavior of the disordered system. According to the theo- 
retical predictions |46I22| , the pseudocritical temperatures 
T| of the disordered system are distributed with a width 
(5[TJ]o„, that scales, in the case of a new random fixed 
point, with the system size as 

5{[T*zU) ^ L-^/r (6) 

In the insets of Figures [4| and [5| we plot these sample- 
to-sample fluctuations of the pseudocritical temperature 
of the magnetic susceptibility. The solid lines show, in 
both cases, a very good power-law fitting giving the values 
0.688(15) and 0.679(21) for the exponent which is also 
in very good agreement with the values denoted above and 
obtained via the classical shift behavior and the accurate 
estimates in the literature |21I25] . 

We move now to investigate the magnetic properties 
of the model. In Figure [HI we provide estimates for the 
magnetic exponent ratio ^ /v of the 3D random-bond Ising 
model, by presenting the FSS behavior of the maxima of 
the disorder-averaged magnetic susceptibility [xlai, for all 
the values of the disorder strength considered. The solid, 
dotted, and dashed lines show a simultaneous fitting of 
the form 

ML - L-^'^ (7) 

using the L = 16 — 64 lattice-range spectrum, providing 
an estimate 1.964(9) for the exponent ratio ^/v- We re- 
mind the reader that the respective fitting attempts of the 
numerical data for each case separately gave the estimates 
-^/v = 1.962(4), 1.967(3), and 1.969(4), for r = 0.75/1.25, 



Fig. 6. FSS of the maxima of the disorder- averaged magnetic 
susceptibility for all the three values of the disorder strength 
considered. The solid, dotted, and dashed hues show a simul- 
taneous fitting attempt of the form [x]Ju ~ L''^" for the range 
L = 16 — 64. The inset shows the corresponding limiting be- 
havior of the ratio i?[x*]ai. defined in the text. 

0.5/1.5, and 0.25/1.75, respectively. All these sets of val- 
ues are in excellent agreement with the accurate estimates 
1.963(5) of reference [H] and 1.97(2) of reference [5S] for 
the ratio ^/v of the site- and bond-diluted models, re- 
spectively. In the corresponding inset of Figure [6| we try 
to further quantify the behavior of the sample-to-sample 
fluctuations of the model by plotting the noise to signal 
ratio Rz - already introduced in Section [2] - as a func- 
tion of the inverse lattice size, for Z = [x*]at)- Clearly, for 
the present model the limiting value of is nonzero, 

indicating, as expected also for marginal disordered sys- 
tems [22]: a strong violation of self-averaging of the mag- 
netic properties of the 3D Ising model with bond disor- 
der, a property that intensifies with increasing disorder 
strength (r 0). We should comment here that analo- 
gous behavior of the ratio -R[x*]„„ has been observed also 
in several other 3D and 2D random models, some of them 
being the bond- [25^ and site-diluted [21 versions of the 
present 3D Ising model and the random-bond versions 
of the square Blume-Capel and triangular Ising [?3] 
model. As a final remark, we stress that the asymptotic 
behavior of these i?-ratios is expected, from theoretical 
renormalization-group arguments, to follow specific power 
law functions and reach limiting values independent of the 
value of the disorder strength 46 . However, such an anal- 
ysis is by far much more demanding, in terms of numer- 
ical data, as it has already been discussed by Berche et 
al. [2S], and, in any case, goes beyond the scope of the 
present work. 

Respectively, in Figure[7|we plot the disorder-averaged 
magnetization at the estimated critical temperatures of 
the disordered model, as a function of the lattice size L. 
The solid, dotted, and dashed lines show a simultaneous 
fitting following the well-known relation 

[MU{T = T,)^L-^/'', (8) 
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Fig. 7. FSS of the disorder-averaged order parameter at the es- 
timated critical temperatures Tdr) for r = 0.75/1.25, 0.5/1.5, 
and 0.25/1.75. The inset illustrates the FSS of the ma:xima of 
the disorder-averaged inverse-temperature derivative of the ab- 
solute order parameter. In both main panel and inset, lattices 
with linear sizes L — 16 — 64 have been used in the fittings. 



giving an estimate of the order of 0.517(5) for the exponent 
ratio Let us note here that, as in the case of the 

magnetic susceptibility in FigureEl the respective separate 
fitting attempts of the numerical data gave comparable to 
the 0.517 estimates in the regime 0.511 — 0.522. Additional 
estimates for this exponent ratio can be obtained from 
the FSS of the derivative of the absolute order parameter 
with respect to inverse temperature d{\M\) /dK. Thus, in 
the corresponding inset of Figure [7] we plot the data for 
d{\M\) /dK averaged over disorder as a function of L, in 
a double logarithmic scale. Their maxima are expected to 
scale with the system size as [48] 



{d{\M\)/dK)* - 



(9) 



As in the main panel, the three lines shown in the inset 
correspond to linear fittings for the three values of the dis- 
order strength, verifying the estimates of the main panel 
of the figure for P/v, being also in very good agreement 
with accurate estimates 0.519(3) [H] and 0.515(5) [53] for 
the site- and bond-diluted models, respectively. 

Overall, the estimates for the critical exponent ly of the 
correlation length and the magnetic exponent ratios j/iy 
and fi/iy given in this Section, and summarized in Table [1] 
are in excellent agreement with accurate estimates in the 
literature for the respective site- 21.. and bond-diluted 
models [53] , reinforcing the scenario of a single distinctive 
universality class in the 3D Ising model with quenched 
uncorrelated disorder, indicating the presence of a new 
random fixed point. 



Table 1. Summary of critical exponents for the pure and dis- 
ordered 3D Ising model. The values of references [31] , [S] , and 
|25| have been used for the pure, site-, and bond-diluted ver- 
sions of the model. 



Disorder Distribution 


V 




fijv 


None (Pure model) 


0.6304(13) 


1.966(3) 


0.517(3) 


Site dilution 


0.6837(53) 


1.963(5) 


0.519(3) 


Bond dilution 


0.68(2) 


1.97(2) 


0.515(5) 


Bond randomness 


0.685(7) 


1.964(9) 


0.517(4) 



the critical behavior of the 3D Ising spin model embedded 
in the simple cubic lattice, by implementing an efficient 
two-stage entropic simulation scheme based on the Wang- 
Landau algorithm. Our numerical approach enabled us to 
simulate large ensembles of disorder realizations of the 
model for very large lattice sizes (up to iV = 64'^) and 
several disorder strengths, and thus, obtain, via finite-size 
scaling techniques, accurate estimates of all the critical 
exponents describing the phase transition of the model. 
Particular interest was paid to the sample-to-sample fluc- 
tuations of the pseudocritical temperatures of the model 
and their scaling behavior that was used as a successful 
alternative approach to criticality. 

The main conclusion of this study can be synopsized 
as follows: the critical behavior of the 3D Ising model with 
quenched uncorrelated disorder is controlled by a new ran- 
dom fixed point, independently of the way randomness is 
implemented in the system. It should be acknowledged at 
this point that, Ballesteros et al. 21i and Berche et al. [IS] 
were the first scientific groups that strongly supported the 
above view for the present model. However, we believe 
that the results brought forward in this paper, combined 
with the existing ones for the site- and bond-diluted mod- 
els, complete the picture, at least for this specific case of 
the 3D disordered Ising model, contributing further to the 
understanding of the concept of universality in disordered 
spin models. 

Closing, we stress here that it would be interesting to 
study further the universality aspects of even more com- 
plicated spin models, in terms of different disorder dis- 
tributions and couplings. One interesting case would be 
to consider field-disorder coupling to the local order pa- 
rameter, and one prominent candidate for this case is the 
3D random-field Ising model [49.50^ . Of course for this 
type of model one may be particularly careful, since in 
this case there exist the well-known hyperscaling violation 
and new scaling concepts should be taken under serious 
account [49 50 51 . Several attempts towards this direc- 
tion of understanding universality in random systems are 
currently under consideration on both numerical and the- 
oretical grounds. 
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